***************************************************************************
***		Do-File to plot exported results and descriptive statistics		***
***		including those computed from SIAB								***
***************************************************************************
clear
set more off
global directory = "INSERT YOUR DIRECTORY"

************************************************************
*** Effect of rint_deuis on profit situation by industry ***
************************************************************

import delimited "$directory\Daten\170204\04_earnsit_lincom.txt", delimiter(space, collapse) varnames(1) clear
gen ind_rank = .
replace ind_rank = 7 if _w73_1 == 1
replace ind_rank = 1 if _w73_1 == 2
replace ind_rank = 2 if _w73_1 == 3
replace ind_rank = 5 if _w73_1 == 4
replace ind_rank = 6 if _w73_1 == 5
replace ind_rank = 4 if _w73_1 == 6
replace ind_rank = 3 if _w73_1 == 7
label define ind_rank 1 "Manufacturing" 2 "Construction" 3 "Services" 4 "Finance" 5 "Retail" 6 "Transportation" 7 "Other"

label value ind_rank ind_rank
gen ind_rank_inv = -1 * ind_rank

qui mean _rint_coef if ind_rank==1
local refline = _b[_rint_coef]

gen estimate = _rint_coef
format estimate %4.3f

twoway  (scatter ind_rank_inv _rint_coef, mlabel(estimate) mlabsize(medsmall) mlabcolor(black) mlabposition(12) msize(vlarge) mlabgap(2) color(black*.7)) ///
		(rcap _rint_cilb _rint_ciub ind_rank_inv, horizontal msize(medlarge) color(black*.7)), ///
/*titles*/			title("") ytitle("") xtitle("") ///
/*legend*/			legend(label(1 "estimate") 	label(2 "95% confidence interval") rows(1) region(lcolor(white)))	///	
/*gen. look*/		plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(black) range(-.5 -7.5)) xscale(lcolor(black) titlegap(2) range(-.125 .075)) ///
					xsize(1) ysize(1.5) ///
/*axis labels*/		xlabel(-.1[.05].05, grid labsize(medsmall)) ///
					ylabel(-1 "Manufacturing" -2 "Construction" -3 "Services" -4 "Finance" -5 "Retail" -6 "Transportation" -7 "Other", notick nogrid angle(0) labgap(*5)) ///
					xline(`refline' , lcolor(black*.3) lwidth(thick)) xline(0, lcolor(black*.05) lwidth(thick))
					graph save Graph "$directory\Output\earnsit_lincom.gph", replace
				

********************************************************************************
***		Prepare SIAB	***
				
use "$directory\Daten\siab_r_7510_v1.dta", clear

* keep only full-time employment spells from Western Germany
keep if quelle_gr == 1 & erwstat_gr==1 & stib!=8 & stib!=9 & stib!=0 & ao_region < 10045

* keep only observations for key date June 30 (doy==181) 
gen begdate = begepi
gen enddate = endepi
format begdate %td
format enddate %td
gen begyear = year(begdate)
gen endyear = year(enddate)
assert begyear==endyear 
gen year = begyear
gen age = year-gebjahr
gen begdoy = doy(begdate)
gen enddoy = doy(enddate)
keep if age >=20 & age<=60 
keep if begdoy <181 & enddoy >=181
bysort persnr year (spell): keep if _n==1		/* for the few pers with multiple spell, keep main spell (see FDZ-Methodenreport 04/2013 for variable spell)  */
bysort persnr year: gen test = _N
assert test==1									/* assert that only one obs per person per year */
drop if year <1993 | year > 2008

* drop public sector 
recode w93_gen_gr (1=1)(2/6=2)(7=7)(8=3)(9=4)(10=5)(11=6)(12=7)(13=7)(14=8), gen(w73_1)
drop if w73_1==8
save "$directory\Daten\temporary.dta", replace

* convert wages into 2010 euros
use "$directory\Daten\deu_cpindex_91-13.dta", clear
merge 1:m year using "$directory\Daten\temporary.dta"
	drop if year <1993 | year > 2008
	erase "$directory\Daten\temporary.dta"
gen rtentgelt = tentgelt_gr * (100/cpindex)
drop if rtentgelt < 20 							/* drop persons who make less than EUR 20 per day, as in LIAB */

*** impute wages
* data prep before impute
gen exper = .
replace exper = age - (6 + 9) 		if bild == 1
replace exper = age - (6 + 9  + 3) 	if bild == 2
replace exper = age - (6 + 13) 		if bild == 3
replace exper = age - (6 + 13 + 3) 	if bild == 4
replace exper = age - (6 + 13 + 5) 	if bild == 5
replace exper = age - (6 + 13 + 5) 	if bild == 6
replace exper = 0 if exper < 0
gen exper2 = exper ^ 2
gen exper3 = exper ^ 3
recode bild (1=1)(2=2)(3=1)(4=3)(5 6=4), gen(bildung) 	/* novtnoabi and novtabi becomes novt, colluni and collfh becomes college */
label define bildung 1 "novt" 2 "vtnoabi" 3 "vtabi" 4 "college"

* federal state
recode ao_region (1000/1999=1)(2000/2999=2)(3000/3999=3)(4000/4999=4)(5000/5999=5)(6000/6999=6)(7000/7999=7)(8000/8999=8)(9000/9999=9)(10000/10999=10), gen(bula)

* identify censored observations
gen highwages = rtentgelt if tentgelt > 80	/* censoring limit is always > 80 */

* censoring limit is the most frequent wage, to account for rounding errors, slightly (3EUR) below --> make sure this also works for east and west berlin!
bysort year: egen bbmg = mode(highwages)
drop highwages
gen cwage = rtentgelt
replace cwage = bbmg - 3 if cwage >= bbmg - 3

* indicator for censored obs
gen cens=.
replace cens = 1 if cwage >= bbmg - 3
replace cens = 0 if cwage < bbmg - 3

* prepare for miimpute command
gen lolim = ln(cwage)
gen uplim = ln(cwage)
replace uplim = . if cens==1

save "$directory\Daten\siab9308.dta", replace

*******************************************************************
***		imputation: separately by education group and industry	***
*******************************************************************
set more off

foreach bild in 1 2 3 4		{
	foreach ind of numlist 1/7		{

		dis " "
		dis "Bildung:" `bild' ", Industry:" `ind'
		dis " "

		* open dataset
		use "$directory\Daten\siab9308.dta" if bildung==`bild' & w73_1==`ind', clear

		mi set wide
		/* impute with force option because very few very small occupations cause perfect collinearity. obs in these occs get miss vals */
		mi impute intreg milnwage i.frau i.deutsch exper exper2 exper3 c.exper#i.frau c.exper2#i.frau c.exper3#i.frau i.bula i.year i.beruf_gr, ll(lolim) ul(uplim) add(1) rseed(239586) force

		replace milnwage = _1_milnwage
		mi unset
		drop milnwage_1_ lolim_1_ uplim_1_ mi_miss lolim uplim
		save $directory\Daten\imp_`bild'_`ind'.dta , replace
		}
	}


***	Stack into one, rename vars, unformat mi ***
set more off
clear

foreach bild in 1 2 3 4		{
	foreach ind of numlist 1/7		{
		qui append using "$directory\Daten\imp_`bild'_`ind'.dta"
		erase "$directory\Daten\imp_`bild'_`ind'.dta"
		}
	}

*** top code at 3x censoring limit ***
replace milnwage = ln(3*(bbmg - 3)) if milnwage > ln(3*(bbmg - 3))

* save
drop _merge
sort year
compress
save "$directory\Daten\siab9308_imputed.dta", replace
erase "$directory\Daten\siab9308.dta"


************************************************************************************************************
***		plot observed and counterfactual profit trends over time for manufacturing and construction	 ***
************************************************************************************************************

import delimited $directory\Daten\04_cftrends.txt, delimiter(space, collapse) varnames(1) clear

twoway 	(line verarb_obs year if year <= 2008, lpattern(solid) lcolor(black*.8) lwidth(thick))  ///
		(line verarb_cf year if year <= 2008, lpattern(shortdash) lcolor(black*.8) lwidth(medthick)), ///
/*titles*/			title("") ytitle("Manufacturing", size(large)) xtitle("") ///
/*legend*/			legend(off)	///	
/*gen. look*/		plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(black) titlegap(5))	xscale(lcolor(black) titlegap(2)) ///
/*axis labels*/		xlabel(none) xtick(1993[1]2008, grid glcolor(white) tstyle(none)) ///
/*grid&ticks*/		ytick(1.25[.25]2.5, grid glcolor(black*.1) glwidth(thin) tlcolor(black) tlwidth(thin)) ///
					ylabel(1.5[.5]2.5, angle(horizontal) grid glcolor(black*.1) glwidth(thin) tlcolor(black) tlwidth(thin) labgap(1)) ///
					ymlabel(1.25 1.75 2.25, angle(horizontal) labcolor(gs4) labsize(small) labgap(2)) ///
					fysize(270) fxsize(120) saving(manufacturing)
		
twoway 	(line bauwi_obs year if year <= 2008, lpattern(solid) lcolor(black*.45) lwidth(thick))  ///
		(line bauwi_cf year if year <= 2008, lpattern(shortdash) lcolor(black*.45) lwidth(medthick)), ///
/*titles*/			title("") ytitle("Construction", size(large)) xtitle("") ///
/*legend*/			legend(label(1 "observed")	label(2 "estimated counterfactual") style(s1mono) rows(1) region(lcolor(white)))	///	
/*gen. look*/		plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(black) titlegap(5))	xscale(lcolor(black) titlegap(2)) ///
/*axis labels*/		xlabel(1995[5]2005, labsize(medsmall)) ///
/*grid&ticks*/		xtick(1993[1]2008, grid glcolor(white) tlcolor(black) tlwidth(thin)) /// 
					ytick(1.25[.25]2.5, grid glcolor(black*.1) glwidth(thin) tlcolor(black) tlwidth(thin)) ///
					ylabel(1.5[.5]2.5, angle(horizontal) grid glcolor(white) glwidth(thin) tlcolor(black) tlwidth(thin) labgap(1)) ///
					ymlabel(1.25 1.75 2.25, angle(horizontal) labcolor(gs4) labsize(small) labgap(2)) ///
					fysize(330) fxsize(120) saving(construction)
		
graph combine manufacturing.gph construction.gph, col(1) ysize(5) xsize(6) imargin(0 0 0 0)	graphregion(color(white))
		graph save Graph "$directory\Output\cfearnsittrends.gph", replace 
		erase manufacturing.gph
		erase construction.gph


********************************************************************************************************
***		plot observed and counterfactual wage trends over time for manufacturing and construction	 ***
********************************************************************************************************

* draw effecton_earn from earnsit_rint and plug into SIAB
import delimited $directory\Daten\170204\04_cftrends.txt, delimiter(space, collapse) varnames(1) clear

sort year
gen effecton_earn_verarb_l0 = verarb_obs - verarb_cf
gen effecton_earn_verarb_l1 = effecton_earn_verarb_l0[_n-1]
gen effecton_earn_verarb_l2 = effecton_earn_verarb_l0[_n-1]
replace effecton_earn_verarb_l1=0 if year==1993
replace effecton_earn_verarb_l2=0 if year==1993 | year==1994

gen effecton_earn_bauwi_l0 = bauwi_obs - bauwi_cf
gen effecton_earn_bauwi_l1 = effecton_earn_bauwi_l0[_n-1]
gen effecton_earn_bauwi_l2 = effecton_earn_bauwi_l0[_n-1]
replace effecton_earn_bauwi_l1=0 if year==1993
replace effecton_earn_bauwi_l2=0 if year==1993 | year==1994

keep year effecton_earn_*
merge 1:m year using "$directory\Daten\siab9308_imputed.dta"
keep if w73_1==2 | w73_1==3						

*** combine coefficients from regwage_reghdfe with effecton_earn to get effecton_lnwage
* manually insert regression coefficients from regwage_reghdfe
global b_ertragslage_l0 = .004628
global b_ertragslage_l1 = .0048434
global b_earnsit_ind_l0 = .0182033
global b_earnsit_ind_l1 = .0142211
global b_earnsit_ind_l2 = .019591
 
gen effecton_lnwage = .
replace effecton_lnwage = effecton_earn_verarb_l0 * $b_ertragslage_l0 ///
						+ effecton_earn_verarb_l1 * $b_ertragslage_l1 ///
						+ effecton_earn_verarb_l0 * $b_earnsit_ind_l0 ///
						+ effecton_earn_verarb_l1 * $b_earnsit_ind_l1 ///
						+ effecton_earn_verarb_l2 * $b_earnsit_ind_l2 if w73_1==2

replace effecton_lnwage = effecton_earn_bauwi_l0 * $b_ertragslage_l0 ///
						+ effecton_earn_bauwi_l1 * $b_ertragslage_l1 ///
						+ effecton_earn_bauwi_l0 * $b_earnsit_ind_l0 ///
						+ effecton_earn_bauwi_l1 * $b_earnsit_ind_l1 ///
						+ effecton_earn_bauwi_l2 * $b_earnsit_ind_l2 if w73_1==3						
						
matrix define cfwagetrends = J(16, 7, .)
matrix colnames cfwagetrends = year verarb_obs verarb_cf N_verarb bauwi_obs bauwi_cf N_bauwi
							
gen cf_lnwage = milnwage - effecton_lnwage
gen mtlmiwage = exp(milnwage) * 30.42 		/* monthly wage */
gen cf_mtlmiwage = exp(cf_lnwage) * 30.42

local yrcount = 0
foreach num of numlist 1993/2008	{
		local ++yrcount
		matrix cfwagetrends [`yrcount', 1] = `num'

		qui mean mtlmiwage if year == `num' & w73_1 == 2	/* manufacturing */
		matrix cfwagetrends [`yrcount', 2] = _b[mtlmiwage]
		qui mean cf_mtlmiwage if year == `num' & w73_1 == 2
		matrix cfwagetrends [`yrcount', 3] = _b[cf_mtlmiwage]
		matrix cfwagetrends [`yrcount', 4] = e(N)

		qui mean mtlmiwage if year == `num' & w73_1 == 3	/* construction */
		matrix cfwagetrends [`yrcount', 5] = _b[mtlmiwage]
		qui mean cf_mtlmiwage if year == `num' & w73_1 == 3
		matrix cfwagetrends [`yrcount', 6] = _b[cf_mtlmiwage]
		matrix cfwagetrends [`yrcount', 7] = e(N)
		}

clear
svmat cfwagetrends, names(col)

twoway 	(line verarb_obs year if year <= 2008, lpattern(solid) lcolor(black*.8) lwidth(thick))  ///
		(line verarb_cf year if year <= 2008, lpattern(shortdash) lcolor(black*.8) lwidth(med)) ///
		(line bauwi_obs year if year <= 2008, lpattern(solid) lcolor(black*.4) lwidth(thick))  ///
		(line bauwi_cf year if year <= 2008, lpattern(shortdash) lcolor(black*.48) lwidth(med)), ///
/*titles*/			title("") ytitle("Euros", size(medium)) xtitle("") ///
/*legend*/			legend(label(1 "manufacturing, observed")	label(2 "estimated counterfactual") ///
						label(3 "construction, observed") label(4 "estimated counterfactual") rows(2) region(lcolor(white)))	///	
/*gen. look*/		plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(black) titlegap(5))	xscale(lcolor(black) titlegap(2)) ///
/*axis labels*/		xlabel(1995[5]2005, labsize(medsmall)) ///
/*grid&ticks*/		xtick(1993[1]2008, grid glcolor(white) glwidth(thin) tlcolor(black) tlwidth(thin)) /// 
					ytick(2800[200]3800, grid glcolor(black*.1) tlcolor(black) tlwidth(thin)) ///
					ylabel(2800[200]3800, angle(horizontal) grid glcolor(white) tlcolor(gs11) tlwidth(thin)) ///
					yline(105, lcolor(white) lwidth(thin))
					graph save Graph "$directory\Output\cfwagetrend_notsplit.gph", replace 


*****************************************************
*** binned scatterplots of effect sizes for 2004 ****
*****************************************************

use "$directory\Daten\siab9308_imputed.dta" if year==2004, replace
save "$directory\Daten\siab04_imputed.dta", replace

* compute effecton_lnwage for all industries
import delimited $directory\Daten\04_earndelta04.txt, delimiter(space, collapse) clear 

gen effecton_lnwage = .		/* compute cumulated effect on wage for the year 2004 for each industry */
foreach num of numlist 1/7 {
				
replace effecton_lnwage = delta04_l0 * $b_ertragslage_l0 ///
						+ delta04_l1 * $b_ertragslage_l1 ///
						+ delta04_l0 * $b_earnsit_ind_l0 ///
						+ delta04_l1 * $b_earnsit_ind_l1 ///
						+ delta04_l2 * $b_earnsit_ind_l2 if w73_1==`num'
	}
keep effecton_lnwage w73_1
merge 1:m w73_1 using "$directory\Daten\siab04_imputed.dta" 	/* plug into 2004 person data */ 
	

* assign persons to the interpercentile range that they fall into based on their avg. real wages
gen milnwage_smear = milnwage + runiform()*.01	/* wage data in SIAB is coarsened, resulting in unequally populated interpercentile ranges. random smearing resolves this */
xtile intpctile = milnwage_smear, nquantiles(100)
bysort intpctile: egen avgeffect = mean(effecton_lnwage)
bysort intpctile: keep if _n==1

*** Effect sizes across 100 interpercentile groups ***
sort intpctile
gen logpoints = 100 * avgeffect
gen logpoints_rel = logpoints - logpoints[50]

scatter logpoints_rel intpctile, msymbol(circle) msize(medlarge) mfcolor(black*.5) mlwidth(thin) mlcolor(black) ///
/*titles*/			title("") ytitle("Log points", size(medlarge)) xtitle("Wage percentile", size(medlarge)) ///
/*legend*/			legend(off)	///	
/*gen. look*/		plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(black) titlegap(5))	xscale(lcolor(black) titlegap(2)) ///
/*axis labels*/		xline(20 50 80, lcolor(black*.1) lwidth(thick)) ///
					xlabel(20 50 80, labsize(medlarge)) /// 
					xmlabel(1 10 30 40 60 70 90 99, notick labsize(medsmall) labgap(2) labcolor(black)) ///
					xline(1 99, lcolor(black*.1) lwidth(thin)) ///
/*grid&ticks*/		xtick(10[10]90, grid glcolor(black*.1) glwidth(thin) tlcolor(black) tlwidth(thin)) /// 
					ylabel(-.6[.2].2, grid glcolor(black*.1) angle(horizontal) tlcolor(black) tlwidth(thin))
					graph save Graph "$directory\Output\percentiles.gph", replace
  

*****************************************************
*** (real log) wage trend by percentile over time ***
*****************************************************
use "$directory\Daten\siab9308_imputed.dta", clear

* matrix with base values only *
matrix define basevals = J(1, 3, .)
_pctile milnwage if year == 1993, p(20 50 80)
	matrix basevals [1, 1] = r(r1)
	matrix basevals [1, 2] = r(r2)
	matrix basevals [1, 3] = r(r3)

* wage growth over time by percentiles
matrix define logtrend = J(16,4,.)
matrix colnames logtrend = p20 p50 p80 year
local yrcount = 0
forvalues year = 1993/2008	{
	local ++yrcount
	dis "yrcount:" `yrcount'
	_pctile milnwage if year == `year', p(20 50 80)
			matrix logtrend [`yrcount', 1] = r(r1) - basevals[1, 1]
			matrix logtrend [`yrcount', 2] = r(r2) - basevals[1, 2]
			matrix logtrend [`yrcount', 3] = r(r3) - basevals[1, 3]
			matrix logtrend [`yrcount', 4] = `year'
		}  
 mat list logtrend 
 clear
 svmat logtrend, names(col)
 
 twoway		 		(line p20 year if year <= 2008, lcolor(black*.4) 	lpattern(solid) 		lwidth(thick)) ///
					(line p50 year if year <= 2008, lcolor(black*.55)	lpattern(longdash) 			lwidth(thick)) ///
					(line p80 year if year <= 2008, lcolor(black*.7) 	lpattern(dash)		lwidth(thick)), ///
/*titles*/			ytitle("") xtitle("") ///
/*legend*/			legend(label(1 "p20") 	label(2 "p50") label(3 "p80") rows(1) region(lcolor(white)))	///	
/*gen. look*/		plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(black) titlegap(5))	xscale(lcolor(black) titlegap(2)) ///
/*axis labels*/		xlabel(1993[2]2008, labsize(medsmall)) ///
/*grid&ticks*/		xtick(1993[1]2008, grid glcolor(white) tlcolor(black) tlwidth(thin)) /// 
					ytick(-.05[.025].1, grid glcolor(black*.1) glwidth(thin) tlcolor(black) tlwidth(thin)) ///
					ylabel(-.05[.05].1, angle(horizontal))
					graph save Graph "$directory\Output\logtrend.gph", replace 

************************************************************************ 
*** Distribution of workers across the wage distribution by industry ***
************************************************************************
use "$directory\Daten\siab9308_imputed.dta" if year==2004, clear
set more off
gen milnwage_smear = milnwage + runiform()*.01	/* wage data in SIAB is coarsened (lumpy), resulting in unequally populated interpercentile ranges. random smearing resolves this */
xtile decile = milnwage_smear, nquantiles(10)


qui tab w73_1 decile, matcell(cells)
	mat list cells

qui tab w73_1, matcell(rowtot)
	mat list rowtot

mat def decilegraph = J(7,10,.)	
forvalues row = 1/7 {
	forvalues col = 1/10 {
		mat decilegraph[`row',`col']=cells[`row',`col']/rowtot[`row',1]
		}
	}
mat list decilegraph

clear	
svmat decilegraph, n(col) 
gen w73_1=_n
reshape long c, i(w73_1) j(decile)
rename c relfreq

graph bar relfreq if w73_1==2, over(decile, label(labsize(large)) axis(lcolor(black*.1))) ///
					title("Manufacturing", size(vlarge) margin(small)) ytitle("")	///
					bar(1, color(black*.4) lcolor(black) lwidth(thin)) ///
					plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(white) titlegap(0)) ///
					ytick(0[.05].2, grid glcolor(black*.1) tlcolor(white) tlwidth(thin)) ///
					ylabel(0[.05].2, labsize(large) angle(horizontal)) yline(.1, lcolor(black*.1) lwidth(thick)) ///
					saving(manufacturing)

graph bar relfreq if w73_1==3, over(decile, label(labsize(large)) axis(lcolor(black*.1))) ///
					title("Construction", size(vlarge) margin(small)) ytitle("")	///
					bar(1, color(black*.4) lcolor(black) lwidth(thin)) ///
					plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(white) titlegap(0)) ///
					ytick(0[.05].2, grid glcolor(black*.1) tlcolor(white) tlwidth(thin)) ///
					ylabel(0[.05].2, labsize(large) angle(horizontal)) yline(.1, lcolor(black*.1) lwidth(thick)) ///
					saving(construction)
					
graph bar relfreq if w73_1==7, over(decile, label(labsize(large)) axis(lcolor(black*.1))) ///
					title("Services", size(vlarge) margin(small)) ytitle("")	///
					bar(1, color(black*.4) lcolor(black) lwidth(thin)) ///
					plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(white) titlegap(0)) ///
					ytick(0[.05].2, grid glcolor(black*.1) tlcolor(white) tlwidth(thin)) ///
					ylabel(0[.05].2, labsize(large) angle(horizontal)) yline(.1, lcolor(black*.1) lwidth(thick)) ///
					saving(services)

graph bar relfreq if w73_1==6, over(decile, label(labsize(large)) axis(lcolor(black*.1))) ///
					title("Finance", size(vlarge) margin(small)) ytitle("")	///
					bar(1, color(black*.4) lcolor(black) lwidth(thin)) ///
					plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(white) titlegap(0)) ///
					ytick(0[.05].2, grid glcolor(black*.1) tlcolor(white) tlwidth(thin)) ///
					ylabel(0[.05].2, labsize(large) angle(horizontal)) yline(.1, lcolor(black*.1) lwidth(thick)) ///
					saving(finance)
			
graph bar relfreq if w73_1==4, over(decile, label(labsize(large)) axis(lcolor(black*.1))) ///
					title("Retail", size(vlarge) margin(small)) ytitle("")	///
					bar(1, color(black*.4) lcolor(black) lwidth(thin)) ///
					plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(white) titlegap(0)) ///
					ytick(0[.05].2, grid glcolor(black*.1) tlcolor(white) tlwidth(thin)) ///
					ylabel(0[.05].2, labsize(large) angle(horizontal)) yline(.1, lcolor(black*.1) lwidth(thick)) ///
					saving(retail)

graph bar relfreq if w73_1==5, over(decile, label(labsize(large)) axis(lcolor(black*.1))) ///
					title("Transportation", size(vlarge) margin(small)) ytitle("")	///
					bar(1, color(black*.4) lcolor(black) lwidth(thin)) ///
					plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(white) titlegap(0)) ///
					ytick(0[.05].2, grid glcolor(black*.1) tlcolor(white) tlwidth(thin)) ///
					ylabel(0[.05].2, labsize(large) angle(horizontal)) yline(.1, lcolor(black*.1) lwidth(thick)) ///
					saving(transportation)				
					
graph bar relfreq if w73_1==1, over(decile, label(labsize(large)) axis(lcolor(black*.1))) ///
					title("Other", size(vlarge) margin(small)) ytitle("")	///
					bar(1, color(black*.4) lcolor(black) lwidth(thin)) ///
					plotregion(color(white))	graphregion(color(white)) ///
					yscale(lcolor(white) titlegap(0)) ///
					ytick(0[.05].2, grid glcolor(black*.1) tlcolor(white) tlwidth(thin)) ///
					ylabel(0[.05].2, labsize(large) angle(horizontal)) yline(.1, lcolor(black*.1) lwidth(thick)) ///
					saving(other)
		
graph combine manufacturing.gph construction.gph services.gph finance.gph retail.gph transportation.gph other.gph, ///
					col(2) ysize(9) xsize(5.5) graphregion(color(white))
graph save Graph "$directory\Output\inddistrib.gph", replace 
	erase manufacturing.gph
	erase construction.gph
	erase services.gph
	erase retail.gph
	erase finance.gph
	erase transportation.gph
	erase other.gph
				
*** END ***
clear
					
					

